$title	ZIP2GDX  -- Read a GTAPAGG zip file and write out data in GDX format

* GTAP7 version. 

$if not set energydata $set energydata yes 

$if not set co2data $set co2data no

$if not set nd $set nd 0

scalar nd	Number of decimals /%nd%/;
abort$(nd<>round(nd)) "Number of decimals must be an integer";


*select dataset
$if not set ds $set ds 10x10



$if not set datadir $set datadir "..\"

*	Extract files from the dataset, retrieve the GDX file
*	and then cleanup:

$if not exist "%datadir%%ds%.zip" $goto missingfile

$if exist "gsdvole.gdx" $call rm gsdvole.gdx

$call 'gmsunzip -C -o "%datadir%%ds%.zip" Basedata.har Default.prm'
$call 'har2gdx Basedata.har'
$call 'har2gdx Default.prm'
$if exist gsdvole.gdx $call 'rm gsdvole.gdx'
$if %energydata%=="yes" $call 'gmsunzip -C -o "%datadir%%ds%.zip"  gsdvole.har'
$if %energydata%=="yes" $call 'har2gdx gsdvole.har'
$call 'gmszip -f -m "%datadir%%ds%" *.har'
$call 'gmszip -f -m "%datadir%%ds%" *.prm'

*up to here ok:

$onecho >setx.gms
$gdxin 'basedata.gdx'
set	i(*)	Goods;
$load i=prod_comm
set x(*)/c,g,i/;
x(i) = yes;
$offecho
$call gams setx gdx=setx

set	x(*)	All goods plus C - G - i and CGDS;
$gdxin 'setx.gdx'
$load x
$gdxin
$call del /q setx.gdx


display x;

*up to here ok.
set	g(x) 	All goods plus C - G - i /c,g,i/,
	i(x)	Market goods,
	f(*)	Factors,
	r(*)	Regions ;




$gdxin 'basedata.gdx'
$load r=reg f=endw_comm i=trad_comm
g(i) =yes;

display g;


*ok


set	src	Sources /domestic, imported/,
	rnum(r)	Numeraire region,
	sf(f)	Sluggish primary factors (sector-specific)
	mf(f)	Mobile primary factors
	ec	Energy goods /
		ecoa	Coal,
		eoil	Crude oil,
		egas	Natural gas,
		ep_c	Refined oil products,
		eely	Electricity
		egdt	Gas distribution /;

alias (r,s), (i,j,jj), (f,ff);

parameters
	vdga(i,r)	Government - domestic purchases at agents' prices,
	viga(i,r)	Government - imports at agents' prices,
	vdgm(i,r)	Government - domestic purchases at market prices,
	vigm(i,r)	Government - imports at market prices,
	vdpa(i,r)	Private households - domestic purchases at agents' prices,
	vipa(i,r)	Private households - imports at agents' prices,
	vdpm(i,r)	Private households - domestic purchases at market prices,
	vipm(i,r)	Private households - imports at market prices
	evoa(f,r)	Endowments - output at agents' prices,
	evfa(f,x,r)	Endowments - firms' purchases at agents' prices,
	vfm(f,x,r)	Endowments - Firms' purchases at market prices,
	vdfa(i,x,r)	Intermediates - firms' domestic purchases at agents' prices,
	vifa(i,x,r)	Intermediates - Firms' imports at agents' prices,
	vdfm(i,x,r)	Intermediates - firms' domestic purchases at market prices,
	vifm(i,x,r)	Intermediates - firms' imports at market prices,
	vims(i,r,s)	Trade - bilateral imports at market prices,
	viws(i,r,s)	Trade - bilateral imports at world prices,
	vxmd(i,r,s)	Trade - bilateral exports at market prices,
	vxwd(i,r,s)	Trade - bilateral exports at world prices,
	vst(i,r)	Trade - exports for international transportation
	vtwr(i,j,r,s)	Trade - Margins for international transportation at world prices
	ftrv(f,j,r)	Taxes - factor employment tax revenue,

	fbep(f,j,r)	Protection - factor-based subsidies,
	isep(i,j,r,src) Protection - intermediate input subsidies,
	osep(i,r)	Protection - ordinary output subsidies,
	adrv(i,r,s)	Protection - anti-dumping duty
	tfrv(i,r,s)	Protection - ordinary import duty
	purv(i,r,s)	Protection - price undertaking export tax equivalent,
	vrrv(i,r,s)	Protection - VER export tax equivalent
	mfrv(i,r,s)	Protection - MFA export tax equivalent
	xtrv(i,r,s)	Protection - ordinary export tax;

$load vdga viga vdgm vigm vdpa vipa vdpm vipm evoa evfa vfm vdfa 
$load vifa vdfm vifm vims viws vxmd vxwd vst vtwr fbep isep 
$load osep adrv tfrv purv vrrv mfrv xtrv 

$load ftrv 

$if errorfree $goto allread

$log *** ftrv not found in HAR file. Value set to zero.
$log *** Ignore the following error message and warning
$clearerror
ftrv(f,j,r) = 0;
$label allread

$gdxin
$call 'rm basedata.gdx'

parameter
	esubd(x)	Elasticity of substitution (M versus D),
	esubva(x)	Elasticity of substitution between factors,
	esubm(x)	Intra-import elasticity of substitution,
	etrae(f)	Elasticity of transformation,
	incpar(x,r)	Expansion parameter in the CDE minimum expenditure function,
	subpar(x,r)	Substitution parameter in CDE minimum expenditure function;

$gdxin 'default.gdx'
$load esubd esubva esubm etrae incpar subpar
$gdxin

$call 'rm default.gdx'


*	Determine which factors are sector-specific 

mf(f) = yes$(etrae(f)=0);
sf(f) = yes$(etrae(f)<0);
display mf,sf;

*	Convert sign of etrae() so that it is non-negative:

etrae(f) = -etrae(f);
etrae(mf) = +inf;

*	Scale the data:

scalar scalefactor /1e-3/;
vdga(i,r) = scalefactor * vdga(i,r);
viga(i,r) = scalefactor * viga(i,r);
vdgm(i,r) = scalefactor * vdgm(i,r);
vigm(i,r) = scalefactor * vigm(i,r);
vdpa(i,r) = scalefactor * vdpa(i,r);
vipa(i,r) = scalefactor * vipa(i,r);
vdpm(i,r) = scalefactor * vdpm(i,r);
vipm(i,r) = scalefactor * vipm(i,r);
evoa(f,r) = scalefactor * evoa(f,r);
evfa(f,j,r) = scalefactor * evfa(f,j,r);
vfm(f,j,r) = scalefactor * vfm(f,j,r);
vdfa(i,x,r) = scalefactor * vdfa(i,x,r);
vifa(i,x,r) = scalefactor * vifa(i,x,r);
vdfm(i,x,r) = scalefactor * vdfm(i,x,r);
vifm(i,x,r) = scalefactor * vifm(i,x,r);
vims(i,r,s) = scalefactor * vims(i,r,s);
viws(i,r,s) = scalefactor * viws(i,r,s);
vxmd(i,r,s) = scalefactor * vxmd(i,r,s);
vxwd(i,r,s) = scalefactor * vxwd(i,r,s);
vst(i,r) = scalefactor * vst(i,r);
vtwr(i,j,r,s) = scalefactor * vtwr(i,j,r,s);
ftrv(f,j,r) = scalefactor * ftrv(f,j,r);
fbep(f,j,r) = scalefactor * fbep(f,j,r);
isep(i,j,r,src) = scalefactor * isep(i,j,r,src);
osep(i,r) = scalefactor * osep(i,r);
adrv(i,r,s) = scalefactor * adrv(i,r,s);
tfrv(i,r,s) = scalefactor * tfrv(i,r,s);
purv(i,r,s) = scalefactor * purv(i,r,s);
vrrv(i,r,s) = scalefactor * vrrv(i,r,s);
mfrv(i,r,s) = scalefactor * mfrv(i,r,s);
xtrv(i,r,s) = scalefactor * xtrv(i,r,s);


*	Drop tiny numbers:

if (nd>0,
	vdga(i,r) = vdga(i,r)$round(vdga(i,r),nd);
	viga(i,r) = viga(i,r)$round(viga(i,r),nd);
	vdgm(i,r) = vdgm(i,r)$round(vdgm(i,r),nd);
	vigm(i,r) = vigm(i,r)$round(vigm(i,r),nd);
	vdpa(i,r) = vdpa(i,r)$round(vdpa(i,r),nd);
	vipa(i,r) = vipa(i,r)$round(vipa(i,r),nd);
	vdpm(i,r) = vdpm(i,r)$round(vdpm(i,r),nd);
	vipm(i,r) = vipm(i,r)$round(vipm(i,r),nd);
	evoa(f,r) = evoa(f,r)$round(evoa(f,r),nd);
	evfa(f,j,r) = evfa(f,j,r)$round(evfa(f,j,r),nd);
	vfm(f,j,r) = vfm(f,j,r)$round(vfm(f,j,r),nd);
	vdfa(i,x,r) = vdfa(i,x,r)$round(vdfa(i,x,r),nd);
	vifa(i,x,r) = vifa(i,x,r)$round(vifa(i,x,r),nd);
	vdfm(i,x,r) = vdfm(i,x,r)$round(vdfm(i,x,r),nd);
	vifm(i,x,r) = vifm(i,x,r)$round(vifm(i,x,r),nd);
	vims(i,r,s) = vims(i,r,s)$round(vims(i,r,s),nd);
	viws(i,r,s) = viws(i,r,s)$round(viws(i,r,s),nd);
	vims(i,r,s) = vims(i,r,s)$viws(i,r,s);
	vxmd(i,r,s) = vxmd(i,r,s)$round(vxmd(i,r,s),nd);
	vxwd(i,r,s) = vxwd(i,r,s)$round(vxwd(i,r,s),nd);
	vxwd(i,r,s) = vxwd(i,r,s)$vxmd(i,r,s);

	vst(i,r) = vst(i,r)$round(vst(i,r),nd);
	vtwr(i,j,r,s) = vtwr(i,j,r,s)$round(vtwr(i,j,r,s),nd);
	ftrv(f,j,r) = ftrv(f,j,r)$round(ftrv(f,j,r),nd);
	fbep(f,j,r) = fbep(f,j,r)$round(fbep(f,j,r),nd);
	isep(i,j,r,src) = isep(i,j,r,src)$round(isep(i,j,r,src),nd);
	osep(i,r) = osep(i,r)$round(osep(i,r),nd);
	adrv(i,r,s) = adrv(i,r,s)$round(adrv(i,r,s),nd);
	tfrv(i,r,s) = tfrv(i,r,s)$round(tfrv(i,r,s),nd);
	purv(i,r,s) = purv(i,r,s)$round(purv(i,r,s),nd);
	vrrv(i,r,s) = vrrv(i,r,s)$round(vrrv(i,r,s),nd);
	mfrv(i,r,s) = mfrv(i,r,s)$round(mfrv(i,r,s),nd);
	xtrv(i,r,s) = xtrv(i,r,s)$round(xtrv(i,r,s),nd);
);


parameter
	evf(ec,i,r)	Volume of input purchases by firms (mtoe)
	evh(ec,r)	Volume of purchases by households (mtoe)
	evt(ec,r,r)	Volume of bilateral trade (mtoe);

$gdxin 'gsdvole.gdx'
$load evf evh evt

$if exist gsdvole.gdx $call 'rm gsdvole.gdx'


*	Declare some intermediate arrays which are required to 
*	evaluate tax rates:

parameter	vdm(i,r)	Aggregate demand for domestic output,
		vom(x,r)	Total supply at market prices,
		voa(x,r)	Output at producer prices;

vdm(i,r) = vdpm(i,r) + vdgm(i,r) + sum{x.local, vdfm(i,x,r)};
vom(i,r) = vdm(i,r) + sum(s, vxmd(i,r,s)) + vst(i,r);
voa(x,r) = sum(f, evfa(f,x,r)) + sum(i, vdfa(i,x,r) + vifa(i,x,r));


vom("i",r) = voa("cgds",r);


parameter
	rto(i,r)	Output (or income) subsidy rates
	rtf(f,x,r)	Primary factor and commodity rates taxes 
	rtpd(i,r)	Private domestic consumption taxes
	rtpi(i,r)	Private import consumption tax rates
	rtgd(i,r)	Government domestic rates
	rtgi(i,r)	Government import tax rates
	rtfd(i,x,r)	Firms domestic tax rates
	rtfi(i,x,r)	Firms' import tax rates
	rtxs(i,r,s)	Export subsidy rates
	rtms(i,r,s)	Import taxes rates;

rto(i,r)$vom(i,r) = 1 - voa(i,r)/vom(i,r);
rtf(f,j,r)$vfm(f,j,r)  = evfa(f,j,r)/vfm(f,j,r) - 1;
rtpd(i,r)$vdpm(i,r)   = vdpa(i,r)/vdpm(i,r) - 1;
rtpi(i,r)$vipm(i,r)   = vipa(i,r)/vipm(i,r) - 1;
rtgd(i,r)$vdgm(i,r)   = vdga(i,r)/vdgm(i,r) - 1;
rtgi(i,r)$vigm(i,r)   = viga(i,r)/vigm(i,r) - 1;
rtfd(i,x,r)$vdfm(i,x,r) = vdfa(i,x,r)/vdfm(i,x,r) - 1;
rtfi(i,x,r)$vifm(i,x,r) = vifa(i,x,r)/vifm(i,x,r) - 1;
rtxs(i,r,s)$vxwd(i,r,s) = 1-vxwd(i,r,s)/vxmd(i,r,s);
rtms(i,r,s)$vims(i,r,s) = vims(i,r,s)/viws(i,r,s) - 1;

parameter	pvxmd(i,s,r)	Import price (power of benchmark tariff)
		pvtwr(i,s,r)	Import price for transport services;

pvxmd(i,s,r) = (1+rtms(i,s,r)) * (1-rtxs(i,s,r));
pvtwr(i,s,r) = 1+rtms(i,s,r);

parameter	
	vtw(j)		Aggregate international transportation services,
	vpm(r)		Aggregate private demand,
	vgm(r)		Aggregate public demand,
	vim(i,r)	Aggregate imports,
	evom(f,r)	Aggregate factor endowment at market prices,
	vb(*)		Current account balance;

vtw(j) = sum(r, vst(j,r));
vpm(r) = sum(i, vdpm(i,r)*(1+rtpd(i,r)) + vipm(i,r)*(1+rtpi(i,r)));
vgm(r) = sum(i, vdgm(i,r)*(1+rtgd(i,r)) + vigm(i,r)*(1+rtgi(i,r)));
vim(i,r) = vipm(i,r) + vigm(i,r) + sum(j, vifm(i,j,r));
evom(f,r) = sum(j, vfm(f,j,r));

parameter	theta(i,r)	Final demand value shares;

theta(i,r) = (vdpm(i,r)*(1+rtpd(i,r))+vipm(i,r)*(1+rtpi(i,r))) /
	sum(j,vdpm(j,r)*(1+rtpd(j,r))+vipm(j,r)*(1+rtpi(j,r)));

parameter	eta(i,r)	Income elasticity of demand,
		epsilon(i,r)	Own-price elasticity of demand;

eta(i,r) = 1-subpar(i,r) - sum(j, theta(j,r)*(1-subpar(j,r))) +
	(1/sum(j,theta(j,r)*incpar(j,r))) * 
	  (incpar(i,r)*subpar(i,r) + sum(j,theta(j,r)*incpar(j,r)*(1-subpar(j,r))));
display eta;

epsilon(i,r)$theta(i,r)
	= (2*(1-subpar(i,r)) 
		- sum(j, theta(j,r)*(1-subpar(j,r))) 
		- (1-subpar(i,r)) / theta(i,r) - eta(i,r)) * theta(i,r);
display epsilon;

*	Define a numeraire region for denominating international
*	transfers:

rnum(r) = yes$(vpm(r)=smax(s,vpm(s)));
display rnum;

*	Translate the dataset to a format which drops vdpm, vipm, rtpd, rtpi, 
*	vdgm, vigm, rtgd, rtgi, vpm and vgm:

vdfm(i,"c",r) = vdpm(i,r);
vifm(i,"c",r) = vipm(i,r);
rtfd(i,"c",r) = rtpd(i,r);
rtfi(i,"c",r) = rtpi(i,r);
vom("c",r) =vpm(r);

vdfm(i,"g",r) = vdgm(i,r);
vifm(i,"g",r) = vigm(i,r);
rtfd(i,"g",r) = rtgd(i,r);
rtfi(i,"g",r) = rtgi(i,r);
vom("g",r) =vgm(r);

*	Transfer investment good CGDS to "i":
*	(Justin Caron 10.2.2010)

vdfm(i,"i",r) = vdfm(i,"CGDS",r);
vifm(i,"i",r) = vifm(i,"CGDS",r);

vdfm(i,"cgds",r) = 0;
vifm(i,"cgds",r) = 0;

rtfd(i,"i",r) = rtfd(i,"cgds",r);
rtfi(i,"i",r) = rtfi(i,"cgds",r);

rtfd(i,"cgds",r) = 0;
rtfi(i,"cgds",r) = 0;
 
vb(r) = vpm(r) + vgm(r) + vom("i",r) 
	- sum(f, evom(f,r))
	- sum(j,  vom(j,r)*rto(j,r))
	- sum(x,  sum(i, vdfm(i,x,r)*rtfd(i,x,r) + vifm(i,x,r)*rtfi(i,x,r)))

	+ sum(i, vdfm(i,"c",r)*rtfd(i,"c",r) + vifm(i,"c",r)*rtfi(i,"c",r))
	+ sum(i, vdfm(i,"g",r)*rtfd(i,"g",r) + vifm(i,"g",r)*rtfi(i,"g",r))

	- sum(j,  sum(f, vfm(f,j,r)*rtf(f,j,r)))
	- sum(i, vdpm(i,r)*rtpd(i,r) + vipm(i,r)*rtpi(i,r))
	- sum(i, vdgm(i,r)*rtgd(i,r) + vigm(i,r)*rtgi(i,r))
	- sum((i,s), rtms(i,s,r) *  (vxmd(i,s,r) * (1-rtxs(i,s,r)) + sum(j,vtwr(j,i,s,r))))
	+ sum((i,s), rtxs(i,r,s) * vxmd(i,r,s));

vb("chksum") = sum(r, vb(r));
display vb;



parameter
	eco2(i,*,r)		Emissions of CO2 (Gg),
	evd(i,*,r)		Volume of energy purchases (mtoe);


eco2(i,g,r) = 0;
evd(i,g,r) =0;

if (nd>0,
	execute_unload '%datadir%%ds%_%nd%.gdx', r, f, g, i, vfm, vdfm, vifm, vxmd, vst, vtwr, 
		rto, rtf, rtfd, rtfi, rtxs, rtms, esubd, esubva, esubm, etrae, eta, epsilon,  evt, evd, eco2;
else
	execute_unload '%datadir%%ds%.gdx',      r, f, g, i, vfm, vdfm, vifm, vxmd, vst, vtwr, 
		rto, rtf, rtfd, rtfi, rtxs, rtms, esubd, esubva, esubm, etrae, eta, epsilon,  evt, evd, eco2;
);



$exit
$label missingfile
$log	*** Error *** Missing file.
$log.
$log	Need to have the following file for this program:
$log.
$log		%datadir%%ds%.zip
$log.
$log	This file is expected be produced by gtapagg.


